c
c     Cosmology and Power Spectrum
c
      subroutine init_csm(cch1,cch2,z) ! __ Set up Cosmological Params
c     cch is a character string of the form h0,Omega0,Omaga_Lambda
      character cch1*(*),cch2*(*)
      real h0,ome0,omel
      common/cospar/om_m,om_l,om_b,h,p_index,sig8
      common/cosmo_today/om_m0,om_l0,sig80
      common/norm/pfac

      if (index(cch1,'- ').le.0) then
         read (cch1,*) h,om_m0,om_l0
         if (abs(om_m0+om_l0-1.).gt.1E-6) then
            om_l0 = 1.-om_m0
            print '(A,A)', '! Only flat universe supported.'
        end if
      else
         h=0.7
         om_m0=0.3
         om_l0=0.7
      endif
c
      if (index(cch2,'- ').le.0) then
         read(cch2,*) sig80,om_b,p_index
      else
         sig80=0.85
         p_index=0.97
         om_b=0.045
      endif
      om_m=om_m0
      om_l=om_l0         
      write (6,'(A,6F6.3)') '! h,OmeM,OmeL,sig8_0,Ome_b,p=',
     &        h,om_m0,om_l0,sig80,om_b,p_index      
      sig8=sig80/grow(z)
      om_m=omega_m(z)
      om_l=omega_v(z)
      pfac=1.0 !! Need rhis. This is used in sigint
      pfac=sig8/sigint(8.)
      pfac=pfac*pfac
      write (6,'(A,4F6.3)') '! OmeM(z),OmeL(z),sig8(z),pfac=',
     &        om_m,om_l,sig8,pfac     
      call init_look(z) ! Initialize Lookup table for Power Spectrum?      
      return
      end
c
c
      real function bias(isam,icen,bmdl,rv_max) ! Total for icen=-1
c
c     integrate over mass function (or, rather, nu function),
c     with a nu-dependent weight, to get linear bias
c
      common/cospar/om_m,om_l,om_b,h,p_index,sig8
      common/cosmo_today/om_m0,om_l0,sig80
      logical first_call,file_exist
      character bmdl*(*),tmpfil*20 ! Bias model
      data first_call/.true./
      first_call = first_call.and.(index(bmdl,'only').ge.1)
      rm_min = 1e8
      rm_max = 5e15
c     ....Critial density 3Ho^2/8\pi G in h^2 solar mass/Mpc^3
      dens_cr = 2.78e11 ! h^2 Solar mass/Mpc^3       
      pi=3.141592654
      nint=249
      rm_min = 1e8  ! Minimum/Maximum halo mass considered
      rm_max = 5e15 
      rf_min=(rm_min*3./(4*pi*dens_cr*om_m0))**(1./3.)
      rnu_min=1.686/sigint(rf_min)
      rf_max=(rm_max*3./(4*pi*dens_cr*om_m0))**(1./3.)
      rnu_max=1.686/sigint(rf_max)
c
c if we risk falling off the bottom of the lookup table,
c just integrate over all nu
c
      dnu=(rnu_max-rnu_min)/float(nint)
      delc=1.686
c
c NB 1.686 arguably needs fixing for Omega ne 1 ???? Not really. Good to high precision....
c
      if (index(bmdl,'smt01').ge.1.or.index(bmdl,'SMT01').ge.1) then
         a = 0.707
         b = 0.5
         c = 0.6
      else if (index(bmdl,'tin05').ge.1.or.
     $        index(bmdl,'TIN05').ge.1) then
         a = 0.707
         b = 0.35
         c = 0.80
      end if
      sum1=0.
      sum2=0.
      if (first_call) then
c     Temporary Output nu vs mass
         i=0
         file_exist=.true.
         do while (file_exist.and.i.le.99)
            write(tmpfil,'(A,I2.2,A)') '_tmp_bM_',i,'.dat'
            inquire (file=tmpfil,exist=file_exist)
            i=i+1
         end do
         open (unit=21,file=tmpfil,status='new')
         write (21,'(A)') '!  nu, bias, Mhalo rvir'
         write (21,'(A)') '!   h   OmeM0 OmeL0  Omeb p_idx sig80 bmdl'
c                            12345 12345 12345 12345 12345 12345 1234
         write (21,'(A,6(1X,F5.2),2A5)') '!',h,om_m0,om_l0,om_b,p_index
     $        ,sig80,bmdl
         write (21,'(A)') '!   Mh     r_vir    bias     nu '
         write (21,'(A)') '!h-1Msun   h-1Mpc'
c                           123456789 1234567 1234567 1234567                         
      end if
c         
      do i=1,nint
         rnu=rnu_min+(i-0.5)*dnu
         rv=rvir(rnu)
         rm=rv**3*800*pi*2.78e11*om_m0/3.
         if (icen.eq.-1) then
            wt=wgt_samp(rm,isam,0)+wgt_samp(rm,isam,1)
         else            
            wt=wgt_samp(rm,isam,icen) ! satellite+center
         end if
         if (index(bmdl,'st99').ge.1.or.index(bmdl,'ST99').ge.1) then
c          Bias Sheth & Tormen         
            bb=1+(rnu*rnu/sqrt(2.) -1 +
     $        0.6/(1.+(rnu*rnu/sqrt(2.))**0.3) )/delc
         else ! SMT01 or TIN05
            bb=1.0+1.0/(sqrt(a)*delc)*(sqrt(a)*(a*rnu*rnu)
     $           +sqrt(a)*b*(a*rnu*rnu)**(1.0-c)-
     $           (a*rnu*rnu)**c/((a*rnu*rnu)**c+b*(1-c)*(1-c/2.0)))    
         end if
         if (first_call) write (21,'(1PE9.2,0PF8.2,F8.2,F8.3)')
     $        rm,rv,bb,rnu
         sum1=sum1+f(rnu)*dnu*wt
         sum2=sum2+bb*f(rnu)*dnu*wt
         if (bb*f(rnu)*dnu*wt/sum2.gt.2e-2/real(nint))
     $        rv_max=rv  ! Set Maximum halo scale fro 1-h term
      enddo
      if (sum1.eq.0) then
         bias = 0
      else
         bias = sum2/sum1
      end if
      if (first_call) then
         close (21)
         write(6,'(A,A)') '! **** Created: ',tmpfil
         first_call = .false.
      end if
      return
      end
c
      subroutine mk_pdat(b_lin2,typ,bmdl)
c                        ^^^^^^ bias*bias
c     y_pdat/1: 1-halo term for sample 2:2-halo term 3:1-halo term for mass.
c
      parameter (npdat=5000)
      real*8 pi,yp1,ypn
      real*8 x_pdat(npdat,2),y_pdat(npdat,2),y2_pdat(npdat,2)
      integer n_pdat(2)
      character typ*(*),bmdl*(*)
      character rline*132,tmpfil*20
      logical file_exist
      common/cospar/om_m,om_l,om_b,h,p_index,sig8
      common/cosmo_today/om_m0,om_l0,sig80
      common/pdat/x_pdat,y_pdat,n_pdat,y2_pdat
c     To maximize numerical stability, we have separate limits and
c     numbers btween one-halo and two halo terms
      data rlgkmin1,rlgkmax1/-1.0,2.5/
      data rlgkmin2,rlgkmax2/-3.0,2.5/     
      pi=3.1415926535
c     Convert Matter PS to non-dim Delta^2(k) and regrid to x-pdat
c     Two halo term
      n_pdat(2)=300
      do i=1,n_pdat(2)
         rk=rlgkmin2+(rlgkmax2-rlgkmin2)*(i-1)/real(n_pdat(2)-1) ! Log Binning
         rk=10.**rk
         x_pdat(i,2)=rk
c        Evaluate power spectrum in halo model
c        (P means dimensionless power Delta^2(k))
         y_pdat(i,2)=p_cdm(rk)*b_lin2 ! 2-halo term
         ! access to p_cdm set up internal variables/musrm t come before the first 
         ! call to powint (needed to calculate 1-h term) 
c         y_pdat(i,1)=powint(rk,typ) ! 1-Halo term for sample
      enddo
c
c     One-halo term, finer linear gridding
c
      n_pdat(1)=1000
      do i=1,n_pdat(1)
         rk=rlgkmin1+(rlgkmax1-rlgkmin1)*(i-1)/real(n_pdat(1)-1)
         rk=10.0**rk
         x_pdat(i,1)=rk
c       Evaluate power spectrum in halo model
c       (P means dimensionless power Delta^2(k))
         y_pdat(i,1)=powint(rk,typ) ! 1-Halo term for sample
      enddo
c      Prepare y2-pdat for spline 
c      do i = 1,2
c         call spline (x_pdat(1,i),y_pdat(1,i),n_pdat(i),yp1,ypn,
c     $        y2_pdat(1,i))
c      end do
      if (index(bmdl,'only').gt.1) then ! Do not calc wprp,xi, dump PS
c     Temporary Output Power Spec File
         i=0
         file_exist=.true.
         do while (file_exist.and.i.le.99)
            write(tmpfil,'(A,I2.2,A)') '_tmp_pk_',i,'.dat'
            inquire (file=tmpfil,exist=file_exist)
            i=i+1
         end do
c
         open (unit=21,file=tmpfil,status='new')
         write (21,'(A)') '!  Dimensionless Power Spec Delta^2(k)'
         write (21,'(A)')
     $        '!   h   OmeM0 OmeL0  Omeb p_idx sig80 prof bmdl'
c                12345 12345 12345 12345 12345 12345 1233 1234
         write (21,'(A,6(1X,F5.2),1X,2A5)')'!', h,om_m0,om_l0,om_b,
     $        p_index,sig80,typ,bmdl
         write (21,'(A)') '!  k[Mpc-1]    Delta^2     ih '
c                            1234567890 12345678901234567890
         do ih=1,2
            do i = 1,n_pdat(ih)
               write (21,'(1PE10.3,1X,E10.3,5X,I1)') x_pdat(i,ih),
     $              y_pdat(i,ih), ih
            end do
         end do
         write (21,'(A)') '!'
         write (6,'(A,A)') '! **** Created: ',tmpfil
         close (21)
      end if
      return
      end
c
      function p_look(rk,ih)
      !
      ! Interpolate/Extrapolate Power-spectrum Lookup table
      !
      parameter (npdat=5000)  
      real*8 x_pdat(npdat,2),y_pdat(npdat,2),y2_pdat(npdat,2)
      integer n_pdat(2)
      real*8 pvalue
      common/pdat/x_pdat,y_pdat,n_pdat,y2_pdat
      !                 
      if(rk.ge.x_pdat(1,ih).and.rk.le.x_pdat(n_pdat(ih),ih)) then
         imin =1
         imax = n_pdat(ih)
c         call splint(x_pdat,y_pdat(1,ih),y2_pdat(1,ih),n_pdat,
c     $        dble(rk),pvalue)
c         p_look=pvalue ! To single precision
          call intra_mono_d (x_pdat(1,ih),dble(rk),imin,imax,ff) ! Two-folding method 
          p_look=y_pdat(imin,ih)*(1.-ff)+y_pdat(imax,ih)*ff     c
          do i=1,n_pdat(ih)
             if(rk.ge.x_pdat(i,ih)) j=i
          end do
          frac=(rk-x_pdat(j,ih))/(x_pdat(j+1,ih)-x_pdat(j,ih))
          p_look=y_pdat(j,ih)*(1.-frac)+y_pdat(j+1,ih)*frac
       else
          if(rk.lt.x_pdat(1,ih)) then
            alp=log(y_pdat(2,ih)/y_pdat(1,ih))/
     $            log(x_pdat(2,ih)/x_pdat(1,ih))
            p_look=y_pdat(1,ih)*(rk/x_pdat(1,ih))**alp
         else
            alp=log(y_pdat(n_pdat(ih),ih)/y_pdat(n_pdat(ih)-1,ih))/
     $           log(x_pdat(n_pdat(ih),ih)/x_pdat(n_pdat(ih)-1,ih))
            p_look=y_pdat(n_pdat(ih),ih)*(rk/x_pdat(n_pdat(ih),ih))**alp
         endif         
      endif
      p_look=p_look*exp(-(rk/200.)**4)      
      return
      end     
c
c    ============ Do not use in the current model ==========
c      subroutine sigvint(sigv,z,wtfun)
c      common/cosmo_today/om_m0,om_l0,sig80
c      external wtfun  ! Weight function
c     
c      pi=3.141592654     
c      !nint=499
c      nint=249
c      dens_cr = 2.78e11 ! h^2 Solar mass/Mpc^3       
c      rf_crit=(rm_crit*3./(4*pi*dens_cr*om_m0))**(1./3.)
c      rnu_crit=1.686/sigint(rf_crit)
c      if (rnu_crit.lt.0.3) rnu_crit=0.
cc     
c     if we risk falling off the bottom of the lookup table,
c     just integrate over all nu
c
c      rnumax=20.
c      dnu=(rnumax-rnu_crit)/float(nint)
c      sigtot=0.
c      tot=0.
c      do i=1,nint         
c         rnu=rnu_crit+(i-0.5)*dnu
c         wt=wtfun(rnu)
c         sig=707.0*sqrt(om_m0)*sqrt(1.+z)*rvir(rnu)
c         sigtot=sigtot+sig*sig*f(rnu)*dnu*wt
c         tot=tot+f(rnu)*dnu*wt
c      enddo
c      
c      sigv=sqrt(sigtot/tot)
c
c      return
c      end
c
c   Power Spectrum Integration
c
      real*8 function xiint(func,r)
      parameter (npdat=5000)
      common/xint/ih,rad
      common/pdat/x_pdat,y_pdat,n_pdat,y2_pdat
      external func
      real*8 x_pdat(npdat,2), y_pdat(npdat,2),y2_pdat(npdat,2)
      integer n_pdat(2)
      real*8 func,ans,ans1,ans2,rkmin,rkmax,dr,err      
!     For dqag
      rad = r  ! Do not remove. Needed to pass r to func.
      dr = dble(r) ! double precision
      i=0
      rkmin = x_pdat(1,ih)
      rkmax = x_pdat(n_pdat(ih),ih)
      call integ_per (func,rkmin,rkmax,dble(r),5.0d-4,ans)
c      call qromb (func,rkmin,rkmax,ans)
c      write (0,*) 'r,ih,ans',r,ih,ans
      if (ans.lt.0) then
         xiint=ans
c         xiint = 0
      else
         xiint=ans
      end if
      return
      end
c
      real*8 function f_wprpint(rk)
      
      implicit real*8 (a-h,o-z)
      real r,p_look,rk4
      common/xint/ih,r
      rk4=rk
      pi=3.14159265439793d0
c      f_xiint=p_look(rk4,ih)*(pi/r/rk/rk)*bessj0(r*rk4) !for w_p(r_p)/rp
      f_wprpint=p_look(rk4,ih)*(pi/rk/rk)*bessj0(r*rk4)  !for w_p(r_p)
      return
      end
c
      real*8 function f_xiint(rk)
      
      implicit real*8 (a-h,o-z)
      real r,p_look,rk4      
      common/xint/ih,r
      rk4=rk
      pi=3.141592654
      ifail=0
      f_xiint=p_look(rk4,ih)*dsin(r*rk)/rk/rk/r  ! for xi(r_p)    
      return
      end
c
      subroutine meff(lgme,galdens,isam,icen) 
!     Effective Mass/Com Number density of gals 
!     In halo model
      real*8 c_mpc,solmass
      real lgme
      common/cospar/om_m,om_l,om_b,h,p_index,sig8
      common/cosmo_today/om_m0,om_l0,sig80
c
      pi=3.141592654
      c_mpc = 3.0856e24
      solmass = 1.989e33
c     ....Critial density 3Ho^2/8\pi G in h^2 solar mass/Mpc^3
      dens_cr = 2.78e11 ! h^2 Solar mass/Mpc^3       
      dens_comv= dens_cr*om_m0
      !nint = 499
      nint=249
      rm_min = 1e8  ! Minimum halo mass considered
      rm_max = 5e15
      rf_min=(rm_min*3./(4*pi*dens_cr*om_m0))**(1./3.)
      rnu_min=1.686/sigint(rf_min)
      rf_max=(rm_max*3./(4*pi*dens_cr*om_m0))**(1./3.)
      rnu_max=1.686/sigint(rf_max)
c     if(rnu_min.lt.0.3) rnu_min=0.
c
c     if we risk falling off the bottom of the lookup table,
c        just integrate over all nu
      rnumax=20.
      dnu=(rnu_max-rnu_min)/float(nint)
      rme=0.
      tot=0.
      dens = 0.
c
      do i=1,nint
         rnu=rnu_min+(i-0.5)*dnu
         rv=rvir(rnu)
         rm=rv**3 * 800*pi*2.78e11*om_m0/3.
         if (icen.eq.-1) then
            wt=wgt_samp(rm,isam,0)+wgt_samp(rm,isam,1)
         else
            wt=wgt_samp(rm,isam,icen)
         end if
         rv=rvir(rnu)
         rmass=rv**3 * 800*pi*dens_comv/3. ! in M_sun
         rme=rme + rmass * f(rnu)*dnu * wt
         tot=tot+f(rnu)*dnu * wt           
      enddo
 10   continue
      if (tot.gt.0) then
         lgme=log10(rme/tot)
         galdens = tot*dens_comv/1e10 !wgt_samp is given weight per 1e10 h^{-1}Msun
      else
         lgme=0.0
         galdens=0.0
      end if
      return
      end
c
c    Subroutines used to calculate Power Spectrum
c    for calculating HOD model. Form J. Peacock's
c     w_aaomega_lrg_org.f...   
c
      function p_cdm(rk)

      common/cospar/om_m,om_l,om_b,h,p_index,sig8

      p_cdm=p_full(rk)

      return
      end


      subroutine init_look(z)

      common/look/rnu_look(1000),rvir_look(1000),con_look(1000)

c
c cover a log grid in filter radius
c

      do i=1,1000

         rfil=-2. + 4.*(i-1)/999.
         rfil=10.**rfil
         rvir_look(i)=rfil/(200.)**(1./3.)
         rnu_look(i)=1.686/sigint(rfil)
c
c       NB 1.686 arguably needs fixing for Omega ne 1
c        No, it is good enough approx
         z_f = z_form(rnu_look(i))
c z_f is evaluated relative to the eopch of observation,
c so need to translate to real z

         z_f = -1 + (1+z_f)*(1+z)
c want delta to be the density contrast relative to the mean at 
c the time of interest    
         dc=3000.* (1.+z_f)**3 / (1.+z)**3
c
c approximate delta => c inversion.
c
c can iterate dc to get consistent halo profile, but wk only cares about c
c
         c=0.282 * dc**0.3333 * ( 1 + 0.5*log10(dc**0.33333 / 36.))
         c=(c/1.7)**0.9
c
c        NFW++
c
         con_look(i)=c

      enddo
      return
      end
c
c
c    Change the input definition on wk from k4,y4 to 
c    rnu,rkv
c
      double precision function wk(rnu,rkv,typ)
c
c     Use NFW (default) or Moore+99, MN310,1147
c
      implicit real*8 (a-h,o-z)
      real rnu, rkv 
      character*3 typ ! 'm99' or 'nfw'
      external f_m99ft, f_nfwft
      common /pro_ft/rkc
      if (typ(1:1).eq.'-'.or.typ.eq.'   ') typ='nfw'
      pi = 3.14159265358999d0
      yc=con(rnu,typ) !Concentration parameter for M99 prof r_v/r_c 
      rkc = rkv/yc ! k measured in 1/r_c unit
      !
      !nint=499
      !nint=249
      !dy=yc/float(nint)     
      sum=0.d0
      if (typ .eq. 'm99' .or. typ .eq. 'M99' ) then ! Moore99 profile
c         do i=0,nint            ! Integrate up to virial radius for FT
c            x=(i+0.5)*dy
c            sum=sum + dy*x*sin(rkc*x)/rkc/x**1.5d0/(1.d0+x**1.5d0)
c     $           *scoeff(i,nint)
c         enddo
         call integ_per (f_m99ft,0.0d0,yc,rkc,1.0d-3,sum)
         wk=sum/((2.d0/3.d0)*log(1.d0+yc**1.5d0))
      else ! NFW profile
c         do i=0,nint            ! Integrate up to virial radius for FTc
c            x=(i+0.5)*dy
c            !taper = exp(-(x/(yc*2))**4.0)
c            sum=sum + dy*x*sin(rkc*x)/rkc/x/(1.d0+x)/(1.0d0+x)
c     $           *scoeff(i,nint)
c         enddo         
         call integ_per (f_nfwft,0.0d0,yc,rkc,1.0d-3,sum)
         wk=sum/(log(1.0d0+yc)-yc/(1.0-yc)) ! Norm from Cooray & Sheth 02 Eq. 76      
      end if
      return
      end

      real*8 function f_m99ft (x) ! integrand/FT of M99 prof
      implicit real*8 (a-h,o-z) ! x is r in units of r_c      
      common /pro_ft/rkc
      if (x.gt.1e-8) then
         f_m99ft = x*dsin(rkc*x)/rkc/x**1.5d0/(1.d0+x**1.5d0)
      else
         f_m99ft = x**0.5d0/(1.d0+x**1.5d0) 
      end if
      return
      end

      real*8 function f_nfwft (x) ! integrand/FT of NFW prof
      implicit real*8 (a-h,o-z) ! x is r in units of r_c      
      common /pro_ft/rkc
c     f_nfw_ft = x*sin(rkc*x)/rkc/x/(1.d0+x)/(1.0d0+x)
      f_nfwft = sin(rkc*x)/rkc/(1.d0+x)/(1.0d0+x)
      return
      end
           
      function f(x)
c     From Sheth, Mo & Tormen 2001 Eq. 6 [f(nu)] OK
c      
      f=0.21617*(1+(sqrt(2.)/x/x)**0.3)*exp(-x*x/2./sqrt(2.))

      return
      end

      function rvir(rnu)

      common/look/rnu_look(1000),rvir_look(1000),con_look(1000)
      imax=1
      do i=1,999
         if(rnu.gt.rnu_look(i)) imax=i
      enddo

      frac = (rnu-rnu_look(imax))/(rnu_look(imax+1)-rnu_look(imax))
      rvir = frac*rvir_look(imax+1)+(1-frac)*rvir_look(imax)

      return
      end
c
c
      function con(rnu,typ)
      character typ*(*)
      common/look/rnu_look(1000),rvir_look(1000),con_look(1000)

      imax=1
      do i=1,999
         if(rnu.gt.rnu_look(i)) imax=i
      enddo

      frac = (rnu-rnu_look(imax)) / (rnu_look(imax+1) - rnu_look(imax))
      con = frac*con_look(imax+1) + (1-frac)*con_look(imax)
      if (typ.eq.'nfw'.or.typ.eq.'NFW') then
         con = con**(1/0.9)*1.7
      end if
      return
      end


      function rn_cdm(rk)
      
      y=p_cdm(rk*0.5)
      yplus=p_cdm(rk*0.5*1.01)
      rn_cdm=-3.+log(yplus/y)*100.5
      
      return
      end
      

      function fnl(y,rn)

      common/cospar/om_m,om_l,om_b,h,p_index,sig8

      g=(5./2.)*om_m/(om_m**(4./7.)-om_l+(1+om_m/2.)*(1+om_l/70.))

      a=0.482*(1.+rn/3.)**(-0.947)
      b=0.226*(1.+rn/3.)**(-1.778)
      alp=3.310*(1.+rn/3.)**(-0.244)
      bet=0.862*(1.+rn/3.)**(-0.287)
      vir=11.55*(1.+rn/3.)**(-0.423)

      fnl=y * ( (1.+ b*y*bet + (a*y)**(alp*bet)) /
     & (1.+ ((a*y)**alp*g*g*g/vir/y**0.5)**bet ) )**(1./bet) 
      
      return
      end

      function grow(z)

      common/cospar/om_m,om_l,om_b,h,p_index,sig8

      onow_m=omega_m(z)
      onow_v=omega_v(z)
      
      grow=gg(onow_m,onow_v)
      grow0=gg(om_m,om_l)
      
      grow=(1+z)*grow0/grow
      
      return
      end
      
      
      function omega_m(z)
      
      common/cospar/om_m,om_l,om_b,h,p_index,sig8
      
      aa=1/(1+z)
      
      omega_t=1.0+(om_m+om_l-1.0)/(1-om_m-om_l+om_l*aa*aa+om_m/aa)
      omega_m=omega_t*om_m/(om_m+om_l*aa*aa*aa)
      return
      
      end


      function omega_v(z)

      common/cospar/om_m,om_l,om_b,h,p_index,sig8

      aa=1/(1+z)

      omega_t=1.0+(om_m+om_l-1.0)/(1-om_m-om_l+om_l*aa*aa+om_m/aa)
      omega_v=omega_t*om_l/(om_l+om_m/aa/aa/aa)
      return
      
      end


      function gg(om_m,om_l)

      gg=(5./2.)*om_m/(om_m**(4./7.)-om_l+(1+om_m/2.)*(1+om_l/70.))

      return
      end



      function z_form(rnu)
c
c knows nothing about era of observation, but "grow"
c works using updated values of omega_m etc. at era of interest
c so the returned z is relative to that time
c
      a=1.
 1    a=a/1.01
      z=-1+1/a
      if(grow(z).gt.(1+1/rnu)) then
         z_form=z
         return
      else
         goto 1
      endif

      end


      function p_full(rk)

      common/cospar/om_m,om_l,om_b,h,p_index,sig8
      common/norm/pfac
      real pfac,rk
c
c approximate COBE normalization allowing for p_index
c
      tilt=p_index-1

      if(om_l.gt.0.) then
         dh=1.94e-5*(om_m**(-0.785-0.05*log(om_m)))*
     $        exp(-0.95*tilt-0.169*tilt*tilt)
         else
         dh=1.95e-5*(om_m**(-0.35-0.19*log(om_m)-0.17*tilt))*
     $           exp(-tilt-0.14*tilt*tilt)
      endif
      p_full=dh*dh*(3000.*rk)**(3+p_index)*tk_eh(rk)*
     $     tk_eh(rk)*pfac
      return
      end


      function tk_eh(yy)
      common/cospar/om_mt,om_lt,om_b,h,p_index,sig8
      common/cosmo_today/om_m,om_l,sig80
c
c the astonishing D.J. Eisenstein & W. Hu fitting formula (ApJ 496 605 [1998])
c remember I use k/h, whereas they use pure k
c
c om_m is the total matter density parameter - i.e. CDM + baryons
c
      rk=yy*h      
      e=exp(1.)
      
      thet=2.728/2.7
      b1=0.313*(om_m*h*h)**(-0.419)*(1+0.607*(om_m*h*h)**0.674)
      b2=0.238*(om_m*h*h)**0.223
      zd=1291*(1+b1*(om_b*h*h)**b2)*(om_m*h*h)**0.251/
     $     (1+0.659*(om_m*h*h)**0.828)
      ze=2.50e4*om_m*h*h/thet**4
      rd=31500*om_b*h*h/thet**4/zd
      re=31500*om_b*h*h/thet**4/ze
      rke=7.46e-2*om_m*h*h/thet**2
      s=(2./3./rke)*sqrt(6./re)*
     $     log((sqrt(1+rd)+sqrt(rd+re))/(1+sqrt(re)))
      rks=1.6*( (om_b*h*h)**0.52 ) * ( (om_m*h*h)**0.73 ) *
     $     (1+(10.4*om_m*h*h)**(-0.95))
      
      q=rk/13.41/rke
      
      y=(1+ze)/(1+zd)
      g=y*(-6*sqrt(1+y)+(2+3*y)*log((sqrt(1+y)+1)/(sqrt(1+y)-1)))
      ab=g*2.07*rke*s/(1+rd)**(0.75)
      
      a1=(46.9*om_m*h*h)**0.670*(1+(32.1*om_m*h*h)**(-0.532))
      a2=(12.0*om_m*h*h)**0.424*(1+(45.0*om_m*h*h)**(-0.582))
      ac=(a1**(-om_b/om_m)) * (a2**(-(om_b/om_m)**3))
      
      b1=0.944/(1+(458*om_m*h*h)**(-0.708))
      b2=(0.395*om_m*h*h)**(-0.0266)
      bc=1/(1+b1*((1-om_b/om_m)**b2-1))
      
      f=1/(1+(rk*s/5.4)**4)
      
      c1=14.2 + 386/(1+69.9*q**1.08)
      c2=14.2/ac + 386/(1+69.9*q**1.08)
      tc=f*log(e+1.8*bc*q)/(log(e+1.8*bc*q)+c1*q*q) +
     &     (1-f)*log(e+1.8*bc*q)/(log(e+1.8*bc*q)+c2*q*q)
      
      bb=0.5+(om_b/om_m) + (3-2*om_b/om_m)*sqrt((17.2*om_m*h*h)**2+1)
      bn=8.41*(om_m*h*h)**0.435
      ss=s/(1+(bn/rk/s)**3)**(1./3.)
      tb=log(e+1.8*q)/(log(e+1.8*q)+c1*q*q)/(1+(rk*s/5.2)**2)
      tb=(tb+ab*exp(-(rk/rks)**1.4)/(1+(bb/rk/s)**3))*sin(rk*ss)/rk/ss
      
      tk_eh=(om_b/om_m)*tb+(1-om_b/om_m)*tc
      
      return
      end
      
      function sigint(r)

      implicit real*8 (a-h,o-z)
      real sigint, r, p_cdm, rk4
      
      !nint=499
      nint=249
      
      sum1=0.d0
      do i=1,nint
         t=(float(i)-0.5)/float(nint)
         y=-1.d0+1.d0/t
         
         rk4=y
         d2=p_cdm(rk4)
         
         x=y*r
         w=(3./x/x/x)*(sin(x)-x*cos(x))
         
         sum1=sum1+w*w*d2/y/t/t
      enddo
      
      sum1=sum1/float(nint)
      sigint=sqrt(sum1)
      
      return
      end

      subroutine interp(c,n)
      real c(500,500),x(500,500)
      do i=1,n
         do j=1,n
            x(i,j)=c(i,j)
         enddo
      enddo
      do i=1,500
         do j=1,500
            fx=(i-1)*(n-1)/499.01
            fy=(j-1)*(n-1)/499.01
            ix=int(fx)
            iy=int(fy)
            fx=fx-ix
            fy=fy-iy
            xa=x(ix+1,iy+1)
            xb=x(ix+2,iy+1)
            xc=x(ix+2,iy+2)
            xd=x(ix+1,iy+2)
            c(i,j)=xa*(1.-fx)*(1.-fy)
     &           +xb*fx*(1.-fy)+xc*fx*fy+xd*(1.-fx)*fy
         enddo
      enddo
      return
      end
      
      function powint(rk,typ)
      ! 
      ! Get power spectrum for 1-halo term
      ! Requires modules to provide weights based on
      ! N(M) wgt1(rnu,icen), wgt2(rnu,icen), wgt_cf(rnu,icen)
      !
      ! icen=1 /one gal @ halo center icen=0, two satellites
      character typ*(*)      
      double precision wk,wind 
c      logical first_run, first_loop     
c      data first_run,first_loop/.true.,.true./
      common/cosmo_today/om_m0,om_l0,sig80
      pi=3.141592654
      nint = 499
      dens_cr =  2.78e11 ! Critical density in h^2 M_sun/Mpc3
      rm_min = 1e8  ! Minimum halo mass considered
      rm_max = 5e15 ! Max mass 
      rf_min=(rm_min*3./(4*pi*dens_cr*om_m0))**(1./3.)
      rnu_min=1.686/sigint(rf_min)
      rf_max=(rm_max*3./(4*pi*dens_cr*om_m0))**(1./3.)
      rnu_max=1.686/sigint(rf_max)
c      if(rnu_min.lt.0.3) rnu_min=0.
c     
      dnu=(rnu_max-rnu_min)/real(nint)
c      
      powint=0.
      tot1=0.
      tot2=0.
      totc = 0.
!
      do i=0,nint         
         rnu=rnu_min+i*dnu
         rv=rvir(rnu)
         rm=rv**3*800*pi*dens_cr*om_m0/3.         
         wt1s=wgt_samp(rm,1,0)
         wt1c=wgt_samp(rm,1,1)
         wt2s=wgt_samp(rm,0,0)
         wt2c=wgt_samp(rm,0,1)
         wt_cf_s = wgt_cf(rm,0) ! Two satellites
         wt_cf_c = wgt_cf(rm,1) ! One center/the other satellite
c
         if (wt_cf_s+wt_cf_c.gt.0) then
            rkv=rk*rvir(rnu)
            win=wk(rnu,rkv,typ)
            powint=powint+f(rnu)*dnu*rkv*rkv*rkv*
     $           (win*win*wt_cf_s+win*wt_cf_c)
            ! (sat-sat)+(sat-center)
            totc=totc+f(rnu)*dnu*(wt_cf_s+wt_cf_c)
         end if
c         if (rk.ge.1.0.and.first_loop.and.rm.ge.1e15) then
c            print *, 'rnu,Mass,powint,rk,win,typ=',rnu,rm,
c     $           powint,rk,win,typ
c            first_loop=.false.
c         end if
         tot1=tot1+f(rnu)*dnu*(wt1s+wt1c)
         tot2=tot2+f(rnu)*dnu*(wt2s+wt2c)
      end do
      first_run = first_loop
      if (tot1*tot2.gt.0) then 
         powint=powint*400./3./3.141592654/tot1/tot2
      else
         powint = 0.0
      end if
      return
      end
c
      real function wgt_samp (rm,isam,icen)
c     Weight applied to each sample, requires galn_halo(rm,isam,icen)
c     icen=0: satellite 
c     icen=1: center    
      rm10 = rm/1e10
      wgt_samp = galn_halo(rm,isam,icen)/rm10      
      return
      end
c     
